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ABSTRACT 

The  nonlinear  Schrodinger  (NLS)  equation  in  a  self-defocusing  Kerr  medium  supports  dark  solitons.  Moreover  the 
mean  field  description  of  a  dilute  Bose-Einstein  condensate  (BEC)  is  described  by  the  Gross-Pitaevskii  equation,  which 
for  a  highly  anisotropic  (cigar-shaped)  magnetic  trap  reduces  to  a  one-dimensional  (ID)  cubic  NLS  in  an  external 
potential.  A  quantum  lattice  algorithm  is  developed  for  the  dark  solitons.  Simulations  are  presented  for  both  black 
(stationary)  solitons  as  well  as  (moving)  dark  solitons.  Collisions  of  dark  solitons  are  compared  with  the  exact  analytic 
solutions  and  coupled  dark-bright  vector  solitons  are  examined.  The  quantum  algorithm  requires  2  qubits  per  scalar  field 
at  each  spatial  node.  The  unitary  collision  operator  quantum  mechanically  entangles  the  on-site  qubits,  and  this 
transitory  entanglement  is  spread  throughout  the  lattice  by  the  streaming  operators.  These  algorithms  are  suitable  for  a 
Type-II  quantum  computers,  with  wave  function  collapse  induced  by  quantum  measurements  required  to  determine  the 
coupling  potentials. 

Keywords:  dark  solitons,  vector  dark-bright  solitons,  nonlinear  Schrodinger  equation,  Gross-Pitaevskii  equation, 
quantum  lattice  representation. 


1.  INTRODUCTION 


The  nonlinear  Schrodinger  equation  (NLS)  plays  a  pivotal  role  in  soliton  physics'  -  whether  for  solitons  propagating  in 
an  optical  fiber  or  for  solitons  in  a  Bose-Einstein  condensate  (BEC).  In  its  simplest  one-dimensional  (ID)  form  the  NLS 
equation  for  the  wave  function  ip  is 

‘1? +  ^ +  oV =  °  (1) 
where  the  potential  F[|  ip  |]  =  1 1 p  |2  for  the  integrable  cubic  NLS.  The  parameter  cr  plays  a  critical  role  in  the  soliton 
properties  :  typically  bright  solitons  exist  for  self-focusing  (cr  >  0)  Kerr  media  while  dark  solitons  for  self-defocusing 
(cr  <  0) .  For  BEC  dynamics  (see  Sec.  4  for  some  more  details),  ip  is  the  macroscopic  wave  function  for  the  condensate 
and  the  two-particle  interaction  strength  a  parameterizes  the  s-wave  scattering  length.  This  scattering  length  describes 
the  elastic  scattering  properties  of  the  atom-atom  interaction  (  o  >  0  for  repulsive  interactions  between  particles  in  the 
condensate,  and  o  <  0  for  attractive  particle  interactions)  and  can  be  readily  tuned  over  orders  of  magnitude  as  well  as 
its  polarity  by  appropriate  variations  in  an  external  magnetic  field  near  the  Feshbach  resonances2  (a  quasi-bound 
molecular  state  at  nearly  zero  energy  coupling  resonantly  to  the  free  state  of  the  colliding  atoms). 


Here,  we  extend  our  earlier  quantum  lattice  algorithms3,4 
possible  solution  to  the  cubic  defocusing  NLS  with  o  v[|  ip  |] 

ip(x,t )  =  axl2Tanf^a}12  [x- xa)^Exj\-2i  at\ 
which  is  a  special  case  of  the  dark  soliton 


ip(x,t)  =  i 


Tanh^a112  b  |x-x0-2^|«|l-b2j  fj 


to  the  case  of  dark  solitons  (i.e.,  o  <  0  ). 
=  -2M  is  the  black  soliton5 


Exp[-2iat] 


The  simplest 

(2) 

(3) 
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with  />-*•!.  The  velocity  of  the  dark  soli  ton  is  .yo(l-/;2j  .  Unlike  the  bright  soliton,  which  has  two  free  parameters, 

the  dark  soliton  only  has  one  free  parameter  that  couples  the  amplitude  to  its  speed  (the  parameter  a  represents  the 
background  amplitude  since  for  dark  solitons,  \ip\  -»  a 1/2  as  |jr|  -*oo.  For  bright  solitons,  \ip\  -*0  as  |  x  |  — *oo). 
The  black  soliton  is  stationary,  with  |  ip  |  -»  0  at  the  center  of  the  soliton  and  a  phase  jump  of  it  across  the  soliton 
minimum 


Fig.  1  A  plot  (both  analytic  and  that  determined  from  our  quantum  algorithm)  of  |  ip  |  for  the  black  soliton  for  the  defocusing  cubic 
NLS.  Black  solitons  are  stationary,  with  zero  intensity  as  its  minimum  and  a  phase  jump  of  it  across  its  minimum.  Dark  solitons 
have  a  non-zero  velocity,  a  non-zero  intensity  minimum  with  a  smaller  phase  jump  (see  Fig.  6). 


Of  more  interest  is  the  collision  of  two  dark  solitons.  An  exact  solution  to  the  cubic  self-defocusing  NLS  Schrodinger 
equation 


i—  +  +  cry[|t/;|]?/>  =  0  ,  with  <7V[|i//|]  =  -2|t//|2 

dt  dx 

that  describes  the  collision  of  two  dark  solitons  of  equal  amplitude  is  given  by5 
tp(x,t)=-  ■ 


a  Cosh  | 

I  +  i  21/2  Sinh\ 

o 

1 

1  Ol 

2UZ  2 mCosh{a2{t- 

-'oil 

j  +  Co.j/z(a 

The  soliton  collision  dynamics  is  shown  in  Fig.  2 


(4) 

(5) 
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Fig.  2  The  time  evolution  of  the  collision  of  two  dark  solitons  of  equal  amplitude  (for  both  analytic  and  the  quantum 
algorithm).  The  spatial  grid  is  1000  and  the  number  of  time  iterations  is  1000. 

Here  we  develop  a  quantum  lattice  algorithm  for  the  evolution  of  dark  solitons  and  their  collisions  as  well  as  the  vector 
NLS  coupling  of  a  dark-bright  soliton  pair  under  a  self-defocusing  nonlinearity.  This  is  of  considerable  interest  since  the 
bright  soliton  only  can  exist  in  a  defocusing  medium  because  it  is  trapped  within  a  co-propagating  dark  soliton.  These 
types  of  interactions  have  been  considered  for  an  inhomogeneous  two-species  Bose-Einstein  condensate.6 


2.  QUANTUM  LATTICE  REPRESENTATION  FOR  NLS 

2.1  Scalar  NLS3’4 

It  is  convenient  to  use  the  one  particle  sector  representation  of  Yepez  and  Boghosian7.  The  spatial  domain  is  partitioned 
into  L  spatial  nodes  :  xe,£  =  l  ....L.  At  each  spatial  site  we  associate  two  qubits,  for  the  scalar  wave  function 

ip  at  xe.  The  simplest  unitary  collision  operator  C  that  will  recover  Eq.  (1)  is  a  factored  operator  that  only  entangles 
the  on-site  qubits  independently  at  each  spatial  site 

C  =  ®i,  Ut  (6) 

A  local  equilibrium  can  be  associated  with  this  on-site  unitary  collision  operator  if  |  v)  is  an  eigenvector  of  U ,  with  unit 
eigenvalue  :  Ut  |  v)  =  |  v) . 
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Introducing  the  -JSWAP  -gate  on  a  site  by  site  basis 
(\  0  0  (A 


0  — -  — -  0 

2  2 

2  2 

0  0  0  1 


with  U,=U 


with  the  Hamiltonian  tensor  product  representation 


U  =  exi 


in  in  t 

1T]exprTV 


a\  o2  +  al  ai  +  o\ 

xx  y  y  z 


where  the  o'  s  are  the  usual  Pauli  spin  matrices  for  the  two  on-site  qubits 
.  (0  l\  .  (0  -i)  .  (I  0) 

oj  ’/H  oj  ’  CT*  "(o  -i]  •  (9) 

It  should  be  noted  that  U4  is  the  identity  matrix,  so  that  U4 1  v)  =  |  v) . 

The  unitary  streaming  operator  St  acts  on  the  qubit  [  at  each  lattice  site  and  shifts  it  to  the  right  to  the  next  spatial 
node.  The  qubit  \q2 )  is  not  streamed  by  Sv  Thus  the  total  collision  matrix  C  (which  independently  couples  the  two  on¬ 
site  qubits  at  each  lattice  site)  does  not  commute  with  the  global  streaming  operators  51,  (which  streams  the  1  st  qubit  at 
each  lattice  site  to  its  nearest  neighbor  on  the  right.  The  algorithm  is  symmetrized  by  also  introducing  the  global 
streaming  operator  S2  which  performs  a  shift  to  the  right  of  the  2nd  qubit. 

The  appropriate  collide-stream  sequence  of  unitary  operators  needed  to  recover  the  dispersive  part  of  NLS  is  dictated  by 
the  fact  that  C4  =  / 

ip(t  +  At)  =  S2  CS2C  S2  C  S2  C.Sf  CS,C  S{  CSyC  ip{t)  (10) 

where  Sj  is  the  transpose  of  5,,  with  Sj  St  =  I  ,  *  =  1,2.  Note  that  the  potential  in  Eq.  (1)  can  be  readily  introduced 
into  the  formalism  by  introducing  a  local  phase  change  to  the  system  wave  function  for  each  on-site  qubit  pair8,7  after 
performing  the  collide-stream  sequence,  Eq.  (10)  : 

tp  -*  A  r]ip  ,  with  a=±  1.  (11) 

To  proceed  to  the  continuum  limit  one  introduces  the  following  scaling: 
the  spatial  shift  between  neighboring  lattice  sites  =  0[e) 

time  advancement  following  stream-collide-phase  sequence,  Eqs.  (10)  and  (11),  A  f =0(e2  j  (12) 

potential  phase  shift  a[|  ^  |]  =  o(e2  J 

In  the  limit  £  -*■  0,  it  can  be  shown  (using  a  symbolic  manipulation  program  like  Mathematica),  that  the  collide-stream- 
phase  sequence,  Eqs.  (10)  and  (11),  under  the  ordering  Eq.  (12)  yields  the  NLS  equation  -  with  arbitrary  potential  V  - 


.dip  9  ip 


+  =  0  +  C>^£2j  ,  with  a=±  1  . 


2.2  Vector  NLS 

To  recover  the  vector  NLS  equations  appropriate  for  the  study  of  polarized  solitons  and  two-species  BEC 
;  dt/h  ,  d  Vi  I  Li.  fU  _  n 


=  ° 


'if +  +  =  0 
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we  now  introduce  at  each  lattice  site  2  qubits  for  and  2  qubits  for  r/r,.  The  collide-stream  sequence  of  unitary 
operators  is  applied  simultaneously  (and  independently)  to  the  qubit  pairs  \cji,cj2)  and  | q.  ,q4 }  since  the  interaction 
between  the  wave  functions  tp,  and  xp2  occurs  only  through  the  potentials  V,;[jr/>l  \\ip2 1]  : 


Vi(f  +  Af)) 

'sl  CS2C  Si  C  S2  C.sl  ci,  csl  cst  c 

ip2(f  +  At)j 

si cs2csl_  CS2C. sj CSt  csl cslc/ 

with 


£xp[  *  ^  il\n>2  Q A  t\  />,  \ 
Exp [  ia2V'2[(t/;1|>|tp2  0A^]j  UJ 


(16) 


3.  NUMERICAL  SIMULATIONS  USING  THE  QUANTUM  LATTICE  ALGORITHM 


3.1  Scalar  Dark  Solitons 

The  quantum  lattice  algorithm  Eqs.  (10)  and  (11),  for  the  case  of  the  black  soliton  solution  to  the  defocusing  NLS 
i*JLt  -  0 

dt  dx2  11 


(17) 


is  shown  in  Fig.  3  for  initial  parameters  a  =  0.01 ,  x0  =  100  and  times  t-\ 00 , 200 .  The  maximal  pointwise  divergence 
from  the  exact  solution,  Eq.  (1),  at  t  =  200  is  typically  less  than  0.05%. 


Black  Soliton  :  I  ip  I 


Fig.  3  Quantum  Lattice  solution  |  T/r  |  of  the  stationary  black  soliton  for  defocusing  NLS.  The  surface  plot  is  shown  in  Fig.  1 
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The  corresponding  quantum  lattice  algorithm,  Eqs.  (10)  and  (1 1).  for  the  dark-dark  soliton  collision  is  shown  in  Fie.  4 
for  initial  parameters:  ,r0  =  500  ,  =  -500  and  a  =  -Jl  /  It). 


Dark-Dark  Soliton  Collision 
Quantum  Lattice  Algorithm 


Dark-Dark  Soliton  Collision 
Quantum  Lattice  Algorithm 


- 1  =  400 

- 1=  500 


i ;  /; 

V  ' 


Dark-Dark  Soliton  Collision 
Quantum  Lattice  Algorithm 


Dark-Dark  Soliton  Collision  :  ERROR  at  t  =  400 
_ Quantum  Lattice  -  Exact  Solution 


-0.0002  I - 1 - 1 _ I _ I _ _ | 

350  400  450  500  550  600  x  650  350  400  450  500  550  600  x  650 

(C)  (d) 

Fig.  4  The  quantum  lattice  algorithm  for  the  dark-dark  soliton  solution  for  the  defocusing  NLS,  Eq.  (17)  with  initial  conditions  : 
x0  =  500  ,  t0 =-500,  a=-\l 2/10 .  (a)  Plot  of  |  ip  |  for  t  =  0  and  t  =  300  as  the  dark  solitons  approach  each  other,  (b)  at  t  =  400  and  t 
=  500  there  is  strong  overlap  followed  by  separation  as  the  dark  solitons  pass  through  each  other,  (c)  times  t  =  700  and  t  =  1000. 

(d)  the  difference  |^(j,r=400)| _ -  \xp(x,t  =  400)  in  the  absolute  values  of  the  wave  function  at  t  =  400.  The  quantum 

simulation  is  on  a  very  coarse  grid  of  only  a  1000  spatial  nodes  and  the  algorithm  takes  less  than  3  s  to  complete  1000  time  iterations 
on  a  MAC  G5  in  double  precision.  The  surface  plot  of  this  dark-dark  soliton  collision  is  shown  in  Fig.  2  The  symmetry  in  the  error 
plot  (d)  indicates  a  very  balanced  algorithm. 
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As  a  final  example  of  the  use  of  our  quantum  lattice  algorithm  we  consider  the  motion  of  a  coupled  dark-bright  soliton  in 
a  vector  dcfocused  NLS  system'’.  For  self-dc focused  NLS,  the  bright  soliton  per  se  is  not  a  solution,  so  that  this  initial 
profile  will  undergo  spatial  dispersion,  Fig.  5(a).  However  such  a  bright  soliton  can  exist  in  self-defocused  NLS  when  it 
is  trapped  within  a  co-propagating  dark  soliton.  The  vector  self-defocusing  NLS  equations  are 


.  dipD  d  Vo 
dt 

i^  + 


dxl 

d2xpB 


dt 


dx 


-  [kaf  +  \M>d  |2  -  p]v»d  =° 

-  [Kl2  +|V'of  -  n\pB  =° 


(18) 


with  exact  dark  soliton  wave  function  ifiD,  and  bright  soliton  wave  functioning  solutions1 


iM*.  0  =  i  -yf/j  Sin  a  +  -tJu  Cos  a. Tank 


K2(l -Tan2a) 


-xn-  ttcTana 


Vb{x,*)  =  ExP 


t  +  i—^KTana 

V2 


■ 

(  x  V 

Seek 

K 

—j=  -  x0  -  tKTana 

It/2  ) 

(19) 


In  the  terminology  of  BEC,  fi  corresponds  to  the  chemical  potential,  N B  is  the  rescaled  number  of  particles  in  the  bright 
soliton  state,  and  the  soliton  inverse  length  is 

Nn 


K  =  J 


2  N  D  IV  D 

a  Cos2 a  +  — —  -  — - 
^  16  4 


The  co-propagating  dark-bright  soliton  pair  has  velocity  -Jl  k  Tana. 


Fig.  5  (a)  For  the  defocusing  scalar  NLS,  the  initial  profile  of  a  bright  soliton  diffuses  since  stationary  solutions  of  the  defocusing 
NLS  are  dark  solitons.  (b)  For  vector  defocusing  NLS,  with  initial  profiles  of  a  bright-dark  soliton  pair  [Eq.  (19)  with  t  =  0],  the 
propagating  dark  soliton  traps  the  bright  soliton  and  prevents  its  decay  even  in  the  defocusing  case.  The  parameters  for  these  runs  are 
/u  =0.000025  ,  N B  =0.001  ,  «=0.5  on  a  spatial  grid  of  14000  and  960  K  time  iterations. 

In  Fig.  6  we  compare  the  real  and  imaginary  parts  of  the  dark  soliton  wave  function  t pD  at  t  =  480  K  and  960  K  to  the 
exact  solution  (19).  There  is  little  error  in  the  Rejt/w  -ipeiact  1  ,  as  seen  Fig.  6  (a).  However,  the  numerical  noise 

L  ^  -1DARK 
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in  Tm  [ipg,  at  (  -  480  K.  Fig.  6  (b),  has  increased  by  two  orders  of  magnitude  into  structured  oscillations 

which  eventually  will  disrupt  the  simulation. 


(a) 


Fig.  6  A  plot  of  (a)  Re[ dark ] ^  and  R^ip dark] exact  at  t  —  480K  and  t  —  960K.  Even  at  t  —  960K,  the  relative  error  is  less  than 

0.07%.  Im [t/»dor*]gL  and  lm[ip  dark]^aa  =  const,  at  (b)t  =  480K,  (c)t  =  960K.  The  extremely  low  numerical  noise  at  t  =  480K 

(errors  <  0.002%)  has  been  amplified  to  structured  oscillations  with  errors  <  0.2%.  The  (theoretical)  phase  discontinuity  of  ipdark 
about  its  minimum  is  0.7 1  Jt  [see  (a)  -  (c)] . 
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4.  DISCUSSION  AND  COMMENTS 


We  have  extended  our  quantum  lattice  algorithm  to  now  cover  both  the  self-focusing  and  self-defocusing  NLS.  While 
this  is  accomplished  by  a  simple  sign  flip  in  the  nonlinear  potential  term,  the  effect  on  the  soliton  solution  is  quite 
dramatic.  In  particular,  the  bright  soliton  complex  wave  function  ipH  (an  exact  solution  to  the  cubic  self-focusing  NLS) 
does  not  have  a  phase  jump  across  its  maximum.  Numerically  this  permits  the  use  of  periodic  boundary  conditions  in 
the  quantum  lattice  algorithm.3'4  On  the  other  hand,  there  is  a  phase  discontinuity  of  tc  at  the  zero  intensity  minimum  of 
a  black  soliton  complex  wave  function.  As  a  result  we  abandoned  the  imposition  of  standard  periodic  boundary 
conditions  and  employed  inflow  streaming  boundary  conditions  using  EOShift 

The  numerical  solution  of  the  defocusing  NLS  for  dark  solitons  is  also  plagued  by  the  occurrence  of  numerical  (round¬ 
off)  instabilities.  Numerical  noise  hits  both  the  standard  split-step  Fourier  method  as  well  as  the  quantum  lattice  method. 
Parameters  have  to  be  chosen  so  that  the  physics  of  interest  is  achieved  before  these  numerical  instabilities  hit. 


The  interplay  between  BEC  and  solitons  has  become  an  extremely  active  research  area.  In  particular,  the  3D  Gross- 
Pitaevskii  (GP)  equation: 


dip(r,t) 


=  -^y2  +  U(r )  +gWc|?p(r,r)|2  ip(r,t) 


is  a  mean  field  theory  of  the  macroscopic  wave  function  of  the  Bose-Einstein  condensate.  U{ r)  is  the  external 

trapping  potential,  g  =  4jrfc2  as  tm  is  the  scattering  amplitude  with  as  the  s-wave  scattering  length,  and  N  is 
proportional  to  the  number  of  condensed  bosons.  The  GP  mean  field  theory  has  been  very  successful  in  describing 
experimental  results  on  dilute  Bose-Einstein  condensates  near  zero  temperature.  Of  much  interest,  both  experimentally 
and  theoretically,  is  the  reduced  dimension  description  of  BEC  as  one  experimentally  constricts  the  condensate  to  disk¬ 
shaped  objects.  In  appropriate  trapping  potentials,  Salasnich  et.  al11  have  derived  a  ID  nonpolynomial  NLS  for  the 
condensate  of  the  form 


h2  d2 


iti^  =  ^  +  V(x) 


271  Jl  +  2asN\(P\ 


r+^L^/1  +  2  asN\<t>\2 


h  +  2asN\<j>\2 


<t>  (21) 


where  ax=-^hlmc is  the  oscillator  strength  in  the  transverse  direction  and  ip(r,t)  =  il)L(y,z,t)<p{x,t).  In  the  special 
case  of  weakly  interacting  Bose  particles,  as  N 1 0 12 « 1 ,  Eq.  (21)  reduces  to  a  cubic  NLS  with  external  potential  V: 
dip  \  h2  d2  Tr/  \  _  sN  .  ,  .2!  , 


ih^-  =  -2—^+  v(x)  +-^r|0f  <P 
dt  2m  dx 2  ;  2jta2'  J 

while  in  the  limit  of  strongly  interacting  axial  high-density  condensate  asN  |0|2»  1  (but  with  overall  3D  dilute 
condensate  as  N  |  ip  |2  «  1 ),  reduces  to  the  quadratically  nonlinear  NLS 


..  dip  %2  d2  \  4  gN"‘  , 

dt  2  mdx2  V  1  2  2nal^2al'  1 

We  intend  to  consider  these  equations  as  well  as  the  full  3D  GP  equation  using  quantum  lattice  algorithms. 


n»JL  =  +  v(x) 
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